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Abstract 

We discuss a new non-linear PDE, u, + (2u xx /u)u x = eu xxx , invariant under scaling of depen- 
dent variable and referred to here as SIdV. It is one of the simplest such translation and space-time 
reflection-symmetric first order advection-dispersion equations. This PDE (with dispersion co- 
efficient unity) was discovered in a genetic programming search for equations sharing the KdV 
solitary wave solution. It provides a bridge between non-linear advection, diffusion and disper- 
sion. Special cases include the mKdV and linear dispersive equations. We identify two conser- 
vation laws, though initial investigations indicate that SIdV does not follow from a polynomial 
Lagrangian of the KdV sort. Nevertheless, it possesses solitary and periodic travelling waves. 
Moreover, numerical simulations reveal recurrence properties usually associated with integrable 
systems. KdV and SIdV are the simplest in an infinite dimensional family of equations sharing 
the KdV solitary wave. SIdV and its generalizations may serve as a testing ground for numerical 
and analytical techniques and be a rich source for further explorations. 

Keywords: genetic programming, advection dispersion equation, travelling waves, recurrence 

1 Introduction 

Several equations of physics have been discovered by searching for one that admits a particular type 
of solution that was known to exist on physical grounds. Perhaps the most celebrated example is 
the discovery of Schrodinger's equation of quantum mechanics via de Broglie's hypothesis that free 
particles are described by plane matter waves. The KdV equation of fluid flow in a canal HI, 

u t + 6uu x + u xxx - (1) 



'Corresponding author. 



was motivated (and in part discovered) by the search for an equation possessing Russell's Q 'wave 
of translation' as a solution. This is the solitary wave 



u(x, = -sech I -—-(x — ct — xrj) I for c > 0, xq e R. (2) 



2 \ 2 

The KdV equation is the simplest conservative 1 -dimensional wave equation with weak advective 
non-linearity and dispersion. Thus it is widely applicable and has been used to model acoustic soli- 
tons in plasmas 0|4]], internal gravity waves [5] in the oceans and even blood pressure pulses [0. 
Its significance was greatly amplified by Zabusky & Kruskal's discovery [3 ] that KdV displays the 
Fermi-Pasta-Ulam recurrent behavior and lack of thermalization Q in a spatially periodic domain 
(see also OEl). This was subsequently attributed to the even more remarkable asymptotic superpo- 
sition principle for scattering of KdV solitary waves and existence of infinitely many local conserved 
quantities in involution [10]. Moreover, KdV is the prototype for an integrable non- linear PDE, its 
initial value problem can be solved by the inverse scattering transform IfTTI . 

We were therefore excited to find that the KdV equation is not the only one with the sech 2 solitary 
wave solution. This serendipitous discovery happened in the course of an investigation undertaken 
by two of the present authors |[T2l to improve the efficiency and accuracy of a genetic programming 
(GP) based method (see [A]) to deduce model equations from a known analytic solution [13]. As a 
benchmark exercise to test out the method for application to nonlinear PDEs, the travelling wave ([2]) 
was given to the program, expecting it to find the KdV equation. But, surprisingly, before finding the 
KdV equation, it found a different equation 



u t + u x = u xxx (3) 



which had the same solitary wave solution. Subsequently, we found that Q is the simplest in a vast 
family of equations sharing the KdV solitary wave. We think of ([3]) as a non-linear wave equation 
for the dispersive advection of the real wave amplitude u. Unlike the KdV equation, where the 
advecting velocity V = 6u is linear, here it is a quotient V = (2u xx /u). As a consequence Q is 
'scale-invariant' under dilation u — > Au of dependent variable. Scale-invariant advective velocities 
have appeared before e.g., the velocity that is invariant under a rescaling of fields, charges and 
currents in a plasma. We refer to ^ as the ( £ = 1 case of the) SIdV equation^ Like KdV, faster SIdV 
solitary waves are narrower, but due to scale invariance, height and speed are generally unrelated, as 
in classical linear wave equations, but unlike KdV and the non-linear Schrodinger equation (NLSE). 

In ^ the dimension- l?/T coefficients of the dispersive and advective terms are both equal to 
one. This is of course very special. More generally we consider the SIdV equation 

(2au xx \ 

u t + \ \u x = eau xxx . (4) 



\ u 

x — > y/ax eliminates the L 3 /T -dimensional constant a, leaving one dimensionless parameter e 
measuring the strength of dispersion relative to advection. Unlike KdV, where they scale differently 
in x, here both scale as L~ 3 . Though we have not found an experimental system modelled by 
it is just as simple and universal among scale-invariant advection-dispersion equations as KdV is 
among all such equations. For generic e, we have identified 2 conserved densities. This is similar to 



'The acronym SIdV highlighting scale-invariance is related to KdV as Sine-Gordon is related to Klein-Gordon. 
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inviscid Eulerian hydrodynamics, but unlike KdV and NLSE which possess an infinite number. While 
SIdV shares solitary waves with KdV at e - 1 , non-linear diffusion and solvability emerge elsewhere. 
Despite being non-linear, scale-invariance ensures SIdV has exact plane wave solutions. Furthermore, 
it possesses bounded spatially periodic travelling waves and similarity solutions. However, SIdV 
cannot arise from a polynomial Lagrangian in the sort of variable that works for KdV. We evade this 
obstruction at some special values of e. Remarkably, numerical evolution of Q shows Fermi-Pasta- 
Ulam-like Birkhoff recurrence @ despite no sign of soliton scattering! 



2 General properties of the SIdV equation 

2.1 Symmetries and Conservation laws 

SIdV ([4]) is a non-linear advection-dispersion equation for the real wave amplitude u(x, t) that is being 
advected by the flow V = 2u xx /u. It is scale-invarian^] under u — > Au. In fact, Q is the simplest 
non-linear translation and scale-invariant advection-dispersion equation that is first order in time. The 
lowest order dispersive term is u xxx , so any such equation can be written as u t + Vu x = aeu xxx for 
some scale-invariant advective velocity V. Without requiring scale-invariance, the simplest choice 

V oc u leads to KdV. KdV is symmetric under space-time (PT) reflection (x, t) — > (-x, -t). Now 
if we also require scale-invariance, the simplest advecting velocities that preserve PT symmetry are 

V oc ^si anc i jl xhe former leads to (pfl) in units where a = 1 . The latter choice too has some notable 
properties (see ^7]). 

Non-zero constants are the simplest solutions of Q. Linearization about a constant yields a plane 
wave with cubic dispersion a> = ek 3 and the characteristic property that the ratio of phase to group 
velocity is 1/3, as in the linear KdV equation. Remarkably, despite being non-linear, Q also admits 
exact plane wave solutions u = A sin (kx - (e - 2)k 3 t + tp] . 

KdV is invariant under Galilean boosts x — > x + ct if u transforms as u — > u(x + ct,t) - (c/6). 
However, SIdV is not invariant under 8u - b+ctu x for any constants b and c t 0. However, Galilean 
boosts could be implemented in a more intricate manner that we have not identified. 

Multiplying by u, and using uu^ x = (uu xx - ^u^j , &4) is written in conservation form 



1 + ^) U X - 6 »»AV. 



0. (5) 



It follows that / = J u 2 dx is conserved. This suggests u 2 is the concentration of some substance 
whose total amount is conserved. Similarly, multiplying Q by u' 2 ^ we get 

(u 1 - 2/£ \ + (2-e)(u- 2/£ u xx ) x = 0. (6) 

So / = f u l ~ 2 ^ dx is also conserved. For e.g., when e = 1, J = J \ dx, so on a bounded domain, J 
is finite for any strictly positive/negative initial condition. These integrals of motion, travelling waves 



(§3.2 1 and numerical evolution (351) indicate that SIdV is generically non-dissipative. 



SIdV &4t is not invariant under dilations of x and has a dimensional scale a which we set to 1 . 
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2.2 Preservation of positivity of u(x) 



At first sight, it appears that u = is a singular point of Q. But u can vanish at points where u x 
or u xx also vanish, provided u xx u x /u is finite, e.g. the above plane wave vanishes at isolated points. 
Among travelling waves, this is generic, near one where u and u xx have common zeros, there is 
another solution with the same property. However, in the numerical and analytical examples studied, 
if u(x, 0) > 0, it remains positive for t > 0. Let us use the second conserved quantity / = J - dx 
to sketch why this is the case for e = 1 on a bounded domain^] Suppose u(x, 0) > and J is finite 
at t — 0. At t\ > 0, let u(x) develop its first zero, this cannot be a first-order zero as u was strictly 
positive. Then J(t\) = oo, contradicting the constancy of 7! So strictly positive initial data cannot 
develop a zero and therefore must remain positivd^ 

2.3 Behaviour of SIdV at some special values of dispersion coefficient 
2.3.1 Reduction to a linear dispersive wave equation when e = -2/3 



For e = -2/3, SIdV may be reduced to a linear dispersive equation for u . If we write uu xx 

i 



{uu x ) x - u 2 - \{u 2 ) x - u 2 x in the conservation law form (|5|> of SIdV, we get 



^d t u 2 + d x 



H)»H(" 2 L 

If e - -2/3, we see that p = u 2 satisfies a linear KdV equation 



0. (7) 



2 

Pt + ^Pxxx = 0. (8) 

So for e = - 1 , p = u 2 is a sort of Cole-Hopf transformation that linearizes the equation. The general 
solution is a linear combination of plane waves 



(9) 



p(x, t) = u 2 (x, t) = f p(k) e ! '( fc *+!* 3 ^ where p(k) = f u 2 (x, 0) e' ikx dx. 
Jr 2/t J r 

This solution could serve as the th order of a perturbative solution of SIdV for nearby e. 
2.3.2 The dispersionless limit e - 

The special case e = of Q gives a non-dispersive non-linear advection equation 

iii + '—^u x - 0. (10) 
u 

It may be written as a conservation law (u 2 ) t + {^ u x ) = f° r me 'charge' density u 2 with flux u 2 . 
Being a second order parabolic PDE, (fl0| may also be viewed as an unusual non-linear diffusion 



3 A similar argument can be given for any e for which a conserved density diverges at u = sufficiently fast. 

4 We assume that u(x) cannot develop an 'integrable' zero for which / is finite. An 'integrable' zero where u ~ \x\" 
for a < 1 would mean u forms a cusp and ceases to be thrice continuously differentiable. Assuming the solution remains 
sufficiently smooth, such possibilities are eliminated. Our numerical simulations did not indicate cusp formation, though it 
is an open question whether SIdV preserves regularity of initial data. 
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equation u t - au xx for the 'temperature' u. The effective thermal diffusivity a = -2u x /u could be of 
either sign. Thus SIdV is a remarkable bridge connecting dispersion, non-linear advection and non- 
linear diffusion. We expect (10 1 to have instabilities if a becomes negative since time-reversed heat 
equations are ill-posed. So we may think of the dispersive term in SIdV as a regularization of ( [TO") ), 
just as KdV is a regularization of the kinematic wave equation (KWE) u t + 6uu x = 0. Remarkably, 
even without a dispersive regularization, (10 1 has smooth solutions. Indeed, unlike the KWE, (10 1 
admits waves that preserve their shape. The general travelling wave of ( 10 1 is u = A cos(kx+2k 3 t+<p) . 
These plane waves may however be unstable, as the diffusivity oscillates in sign! 



2.3.3 SIdV to mKdV when e = 2/3 

We did not find any value of e at which a transformation reduces the SIdV to the KdV equation. But 
there are special values of e at which it comes close. First, the rational non-linearity of SIdV can be 
written as a polynomial while retaining the advection-dispersion structure of the equation. To do so, 
we use the invariance of SIdV under rescaling to choose u dimensionless and put u - e w . Therj^] 

w, + [(2 - e)w 2 x + (2 - 3e)w xx \ w x = ew xxx . (1 1) 

This polynomial form of SIdV was convenient for numerical evolution and also in our search for a 
Lagrangian. It also indicates that e = 2, 1 are somewhat special. At these values, we get KdV-like 
dispersive wave equations with advecting velocities oc w 2 and w xx . These are among the simplest PT 
symmetric advecting velocities beyond KdV. Moreover, at e = 2/3, the sign of the 'local diffusivity' 
is reversed. We see qualitative effects of this reversal in the stability of our numerical simulations as 
e is decreased below 2/3. What is more, at e = 2/3 SIdV is reducible to the modified KdV (mKdV) 
equation. Differentiating in x at e - 2/3, putting v = w x and letting x — > -x and t — » \t, we get the 
defocusing mKdV equation, which is integrable and related to KdV via the Miura transform Ifl4l 

v t - 6v 2 v x + v 3x = 0. (12) 



3 Some similarity and travelling wave solutions of SIdV 
3.1 Similarity solutions 

SIdV is invariant undeij^] x — > Ax and t ^ A 3 t. So we seek solutions u(x, t) - f{z) in the similarity 
variable z - x 3 /54t . We get a third order non-linear ODE for f(z) 

^f'f ~ z 2 f"f + ezf'f - jf' 2 + (| + z)f'f = 0. (13) 
By the substitution g = f If we reduce this to a cubic 2 nd order ODE with variable coefficients 

€ -Y g " + (I ~ l Y gg ' + ezg ' + [l - l ) z2g3 + ( £ - l) zg2 + {l +z ) g = - (14) 

5 Of course, if we restrict to real w, this equation will apply to solutions where u remains everywhere positive. 
6 More generally, SIdV is invariant under x — > Ax, t — > A 3 t, u — > A r u for any y. So r rl3 u(x, t) - f(z) is a scale invariant 
combination. Since y is arbitrary, we restrict here to the simplest case y = . 
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We haven't solved the similarity ODE in general, but in the dispersionless limit e = 0, it becomes a 
linear ODE with a regular singularity at z = 

zf" + \f'-f = 0. (15) 

The two linearly independent solutions may be expressed in terms of the confluent hypergeometric 
function oF\(a,z) or the modified Bessel function of the 1 st kind I n (z): 

Mz)= oFil^zj = 1 + | + = (2^) and 

/ 2 ( Z ) = Z 5 F! I ~,zj - zMl + j + ...J = rl~U/i(2V^). (16) 

Both f\ and /2 are monotonic and grow « e 2 Vz as z — > °°. So they are bounded at late times (z — > 
or f » i 3 ) but are unbounded at early times (z — > °o or jc 3 » ? — > 0). Interestingly, there is a 
unique (up to scale) linear combination f(z) = r(4/3)/i(z)-r(2/3)/2(z) that is bounded for all z > 
(x, t > 0). It begins at f(0) = T(4/3) and monotonically decays to zero as z — > °°. 



3.2 Travelling waves 



Here we discuss travelling waves (u = f{%) with £, = x - ct) for SIdV Q on the unbounded domain 
-oo < £ < oo . Travelling waves must satisfy the third order non-linear ODE -off + Iff" - ef f" = 
0. We may write this in 'conservation law' form 



and integrate once to get 



2eff" -(e + 2)/' 2 + c/ 2 + 3fi = 0. 
The substitutions p - f and F — p 2 give us a first order linear ODE for F{f) : 



(17) 
(18) 



efF'(f) -{e + 2)F + cf + W = 



or 



eF'(w) - (e + 2)F + c<? 2w + 3fi = 0, (19) 



where / = e w . This inhomogeneous 1 st order ODE is reduced to quadrature using the integrating 
factor e - {l+2/£)w . For e * 0, in terms of r(w) = e~ (l+2/£)w F(w) , we get 



er' + 3Be- {1+2/£)w + ce il - 2/£)w = 0. 



For £ ^ fj, ±2 , we integrate to get 



F(w) = Ae( l+ ~> + 



ce e 2w 3B 
+ 



2-6 e+2 

In other words, the equation for travelling waves has been reduced to quadrature 



(20) 



(21) 



df 



where F = 



2 



|A/ 2 -c/ 2 log/+f 



if £ - 0. 
if £ = 2, 



f/ 2 + f logZ-iA if £--2, and 
A/ ( 1 + i) + ( ^ ) / 2 + ^ otherwise. 



(22) 



6 



This travelling wave integral can be understood by a mechanical analogy lfT4l . It is the zero 'energy' 
condition E = f' 2 + V(f) = for the 'coordinate' /(ft at 'time' £ of a non-relativistic particle of 
mass 2 moving in the 1 -dimensional potential V(f) = -F(f) . Bounded travelling waves correspond 
to bound trajectories of this particle. Plotting V(f) shows that, for appropriate ranges of A, B and 
c , there are bounded spatially periodic/solitary travelling waves with heights between successive real 
simple/double zeros of V(f). For generic e, the travelling wave integral J F^ 1 ^ 2 df cannot be eval- 
uated using elementary functions. But for e = oo, 0,-2/3, the integral is trigonometric/exponential 
and for e = 2/3, 1 it is elliptic. 

To illustrate, we consider the case e = 1 where SIdV shares solitary wave solutions with KdV. 
Here we must evaluate an elliptic integral J(Af 3 + cf 2 + B)~ l ^ 2 df, and bounded travelling waves are 
(limits of) cnoidal waves. Suppose A > and F(f) = Af 3 + cf 2 + B = A(f - f)(f - f 2 )(f - fa) has 
three simple real zeros > f\ < f 2 < fj, > . Then we have a periodic cnoidal wave with trough at 
fi and crest at f 2 , determined by 



±(£-£l)= f Ug where /(ft) = /i. (23) 

Jfi 



dg_ 

y/A(g - - f 2 )(g - / 3 ) 



Transforming to g = f\ + (J2 - fi) sin 2 6 and defining the shape parameter < m = y~t < 1 , (|23|> 



/3-/1 



becomes a standard incomplete elliptic integral of the first kind. It is inverted in terms of a Jacobi 
elliptic function cn(w;ra) - cos0, with modulus k (m = k 2 ) 



2u r d6 

Z = ti±^^^= where u — (24) 
y/Mfs-fi) Jo Vl-msin 2 ^ 

At the upper limit g = f = f 2 - (J2 - fi) cos 2 cp, so the cnoidal wave for A > is 



f = h-(fi- fx) cn 2 |i yfMf^fOtf - ft) ; mj . (25) 



Its shape depends on m while its wavelength and speed are 



4K(m) 

A= ; : and c = -A(f 1+ f 2 +f 3 ). (26) 

^Mft-fi) 



Here K(m) = f^ 2 , de is the complete elliptic integral of the 1 st kind. 

JU Vl-'«sin 2 

The advecting velocity field V = -y- for cnoidal waves is finite since / and /" have common 
zeros. Moreover, by modifying the parameters A, ft, we get nearby waves with the same feature. 
Unlike for KdV, where the shape, speed and wavelength of cnoidal waves are non-trivially modified 
upon a rescaling of amplitude, here / h-> Af produces a new cnoidal wave with the same m , c , A and 
phase ft , since the constants transform as (A, B,f) \-> (A/ A, A 2 B, Af). 

If A < , the cnoidal wave extends between f% and fi and is given by 

/(£> = h + (/3 - fi) cn 2 I X - VA(/i-/ 3 ) (f - ft) ; (27) 
where m = (/ 3 - / 2 )/(/ 3 - fx) and A = 4K(m) ^/A(f-f 3 ). 
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Solitary waves are cnoidal waves of infinite wavelength. They occur when a pair of simple zeros 



of F coalesce to form a double zero. For example, if — > f% in (25 1 holding and A > fixed, 



then m — > 1 , A'(m) — > oo and we get a left-moving solitary wave of depression 

M) -> /2 - (72 - /i) sech 2 1 1 V^(/2 - - ft)) for A>0. 



If the zeros f\ -> /2 coalesce in ( 27 1 we get a left-moving solitary wave of elevation 



/(£) -» /2 + (/3 - /2) sech 2 1 ~ V^(/2 - / 3 Xf - for A<0. 



(28) 



(29) 



Finally, if A — > , cnoidal waves reduce to sinusoidal waves with cubic dispersion, 

/(£) - sin ( V 3 ? (M - ft)) . c < 0, arbitrary. (30) 

4 Search for a variational principle for SIdV 

To relate symmetries to integrals of motion, it is interesting to find a Lagrangian or Hamiltonian 
formulation for SIdV. The existence of a Lagrangian for a given equation depends on the field 
variables used and the sort of Lagrangian allowed. For instance, the dispersive term u xxx in KdV 
u, + 6uu x + u xxx = isn't the variation of any polynomial in u and its derivatives. This is because 
every quadratic differential polynomial in u involving three ;t -derivatives is a total derivative. Yet, as 
is well-known, if we put u = Xx > KdV follows from the Lagrangian density \xtXx - x\ ~ \x 2 xx ■ 

In looking for a Lagrangian for SIdV, we choose to work with w = logw, which satisfies the 



KdV-like equation (111 with polynomial non-linearity. By analogy with KdV, we put w = (f> x and 



seek a polynomial action in <p and its derivatives, whose Euler-Lagrange (EL) equations are 

<p xt + [(2 - e)(f> 2 xx + (2 - 3e)0„ x ] <p xx = e<p xxxx . (31) 

<p is a natural variable since the linear part of the equation (f> xt = e<p4 x admits the polynomial La- 
grangian Xo = \<Pt<t>x + § #u- We wish to add potentials V, W to Xo to reproduce the quadratic 



4>xx4>3x and cubic terms (p xx in ( [31b . To do so, we note a couple of general features. If V is a mono- 
mial of degree n in <p and its x -derivatives, then (1) the resulting terms in the EL equation form a 
differential polynomial of degree n — 1, and (2) the total number of ^-derivatives in each term of 
the differential polynomial are the same as in V. Therefore, to produce <p xx in the EL equation, W 
must be a quartic differential polynomial with 6 x -derivatives, and to give <p xx (f>3 x , V must be a cubic 
differential polynomial with 5 derivatives. In other words, 

V = a i <p<p<p5 x + Cl2 (f>(f>x<f>4x + «3 002*03* + «4 0.v0jc03x + «5 0*02*02* = ^ O-Ni- (32) 

i 

Can a\ be chosen so that the variation of J Vdx gives (2 - 3e)0„03^? Unfortunately not, a^] 
S C 

~~r V((f>(y)) dy = - (10ai - 5a 2 + «3 + 4a 4 - 2a 5 ) (2(p xx <f> xxx + X 4 *) . (33) 
x) J 



6<f>( 



7 In particular, this means all the V, differ from one another by total x -derivatives. 
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For no choice of a, can we produce just a quadratic monomial oc <p xx <p xxx in the EL equation. Simi- 
larly, we showed that there is no quartic differential polynomial W that gives <p\ x upon variation. We 
conclude that there is no polynomial Lagrangian in and its derivatives leading to SIdV. 

However, the polynomiality assumption is quite strong. There may be a non-polynomial La- 
grangian in (j> or one in a variable non- locally related to (p . An interesting example of such a pos- 
sibility occurs when e = -| and SIdV becomes the linear dispersive equation p t + (2/3)p xxx = 
upon substituting p = u 2 . This equation follows from X = \tpttyx - jiAL' where if/ x - p. However, 
\p = J x e 2 ^ y dy is non- locally related to (/> , so the Lagrangian is non-local in cp . 

Another way around this negative result is that there may be a Hamiltonian that is a differential 
polynomial in cp, but with non-canonical Poisson brackets. Such a possibility is realized if e = 2/3, 
when SIdV can be transformed into the mKdV equation v t - 6v 2 v x + v^ x = by the substitution v = ^ 
and a rescaling (§ 2.3.3 1. mKdV admits a Hamiltonian formulation v t = {H, v] with 

H= \J( v2x+ j) dX < v «,v(y)} = -d x 5(x-30. (34) 

So for e - 2/3, SIdV admits a polynomial Hamiltonian in the variable v = (logw) v - w x = <p xx . It 
would be interesting to find a Hamiltonian/Lagrangian formulation of SIdV for other values of e. 



5 Numerical evolution of SIdV solitary waves: recurrent behavior 

We numerically solved^Jfhe SIdV initial value problem for e = 1 on the interval [-n, n] with periodic 
boundary conditions. Numerical evolution of one sech 2 wave produced a right-moving travelling 
wave, as expected from the exact solution on (— co, oo). We also considered two solitary wave^] 

When the two solitary waves were initially separated by some distance but had the same heights and 
speed, they were observed to travel without much interaction, just like individual solitary waves. 

Next, we gave the two waves the same initial amplitudes A 12 = 1 but different speeds (or widths) 



(c\,C2) = (4,2) and locations {x\,x-i) = (-n/2,n/4) at t = (fig. 1(a) 1. So it would take each 
wave by itself a time of | and n to traverse the 27r -interval. Solitary waves that decay at 00 must be 
right-moving, so we couldn't give them opposing velocities. But c\ > C2, so wave- 1 caught up with 
wave- 2 due to periodic boundary conditions and collided with it from the rear. Then they separated 
into a small leading wave and a larger trailing wave. The original solitary waves did not retain their 
shapes. The smaller wave that emerged from the collision moved faster and caught up with the bigger 
one by going round the circle. During the next collision, the smaller wave rear-ended the larger one. 
The large one in front morphed into a fast small wave, leaving behind a slow large wave. This is 
illustrated in the last three plots of fig. l(a)| Qualitatively, this pattern seemed to repeat as the IVP 



was solved up to t = 40 , allowing more than two dozen collisions to be observe' 



8 The evolution was done using NDSolve on Mathematica. Stability of the numerical evolution was slightly enhanced 
by working with w = log u which satisfies w, = W3 X + w x w xx — w\ . Positivity of u(x, 0) was preserved at all times. 

'The sum over jeZ ensures that the initial condition satisfies periodic boundary conditions u(-n) = u(n) . In practice, 
the sum was restricted to \j\ < 3 since the remaining terms are exponentially small for —n < x < n. 

10 A maximum of 1200 grid points were placed at an average spacing of 0.1% the domain width (2n). The qualitative 
features reported here were unchanged by adding 200 more points. Over the time interval < / < 40, 1800 nearly equal 
steps were taken with an average time step of 1 .4% of the time it took wave- 1 in isolation to traverse the domain. 
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12 3 




12 3 




2n J-x 



tlx Versus time 



(a) u(x) at t = 0, before, during and after collision of (b) f m 2 dx is conserved upto numerical errors, which 
right-moving waves on 2n -periodic interval. cause it to grow by about .3% for < t < 40. 



Figure 1: SIdV evolution of two solitary wave initial state. Times < t < 40 are indicated above the 



plots in fig. 1(a) 



Though not KdV solitons, the solitary waves that were involved in these collisions displayed a 
certain coherence, they did not dissipate nor degenerate into ripples. Despite not being periodic, the 
evolution seemed to approximately revisit earlier configurations. Interestingly, there was no equipar- 



titioning of wave intensity. We illustrate this in fig. 2(a) by plotting the absolute squares of the first 
few Fourier coefficients of u{x) as a function of time) 11 ! There is some exchange of intensity among 
the first 3 or so Fourier modes £0,1,2, with an approximate periodicity of T « 3. But there is no 
appreciable leakage to higher Fourier modes |c3 4,5| 2 , which are uniformly three to six orders of mag- 
nitude smaller than |co| 2 . There is some growth in |c4,5| . But this can't be distinguished from an 
accumulation of numerical errors, which also caused the integral of motion / = J u 2 dx to increase 
by 0.3% over a time < t < 40 (fig. [T(b)] >. 

What is more, though the Raleigh quotient ('gradient energy' or mean square mode number) 



G(0 



/ !"*(*. Ol 2 d * = z„k 2 |c„(0| 2 
J\u(x,t)\ 2 dx ZJc„(0l 2 



(36) 



is not conserved, it seems to oscillate between bounded limits (fig 2(b) I. The rms mode number 
hovers around v - V<2 ~ \ ■ A calibration of Q using the linear dispersive equation u t - u^ x 
indicates that there are about 2N + 1 = V12<2 + 1 ~ 2 active degrees of freedom present, cq being 
the dominant one, with some contribution from c±\ . These numerical simulations indicate that in a 
periodic domain, the SIdV equation displays recurrent behaviour (for e = 1 ) despite possessing only 
two (known) constants of motion. 



6 Advection-dispersion equations sharing KdV solitary waves 

One of the questions that puzzled us after the discovery of ([3]> by genetic programming, was whether 
there are other such non-linear advection-dispersion equations, sharing the sech 2 wave with KdV 
lfT6l[T71[T8l . To explore this question let us consider the following generalized form of an advection 

n c„(f)= J u(x, t) e~"' x dx, the negative coefficients |c_„| 2 = |c„| 2 contain no new information for real u(x,i). 
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(a) log 10 |c„(f)l 2 for modes < n < 5 (top to bottom). (b) Rayleigh quotient from 5 th order Fourier series 



Figure 2: Time dependence of Fourier coefficients of u(x) and gradient energy for a two-solitary- 
wave initial state. 



dispersion equation, 

iif + Vu x = 8 a u xxx , (37) 

where 8 is an arbitrary dimensionless parameter and V(u, u x . . .) is an arbitrary function. We work 
in units where a = 1 and assume the advecting velocity V(u, u x . . .) to be translation-invariant, so 



that constants are generically solutions. We look for all V and 8 for which (37 ) admits every asymp- 
totically decaying KdV travelling wave as a solution. We proceed by supposing that every decaying 
travelling wave u(x, t) = f(x - ct) = that solves ([T]) also solves ( [37] ) for the same speed c . Then 

- cf + 6ff + f" - and - cf + Vf - 8f" = 0. (38) 

Eliminating /"' we get 

{V-(8+\)c + 68f}f = 0. (39) 

This equation for V and 5 must be satisfied for arbitrary decaying KdV travelling waves. Since KdV 
admits non-constant travelling waves, /' £ 0. So we must have V = (5 + l)c - 66f. c is eliminated 
using the KdV equation -cf + 6ff + f" = 0. However, simply substituting c = 6f + f" If gives 
V\ - 6u + (5 + l)^ 1 leading back to the KdV equation. On the other hand, we could integrate either 
once or twice while omitting integration constants for asymptotically decaying waves and find 

-c/ + 3/ 2 +r = or - C -f+f+ l -(ff = 0. (40) 



Eliminating c from (39), we find 5 is arbitrary and V can take the functional forms 

2 

V z = 3(1 -S)u + {\ + 8)— or V 3 = 2(1 - 28)u + (1 + 8y±. (41) 

u u 1 

However, these are not the only possibilities. Instead of integrating -cf + 6ff + f" - 0, we could 
differentiate this expression (any number of times!) and then use it to eliminate c. In this manner we 
get an infinite sequence of higher order advecting velocities for which ( 37 ) admits ([2]) as a solution: 



Va = 6m + (1+5) ( — £ + — , V5 = 6u + (1 + 8)1 + — 
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(2Aut, x u x 18 

\ M4v U^[ x U4 X ) 



V 6 = 6u + (l + S) p^_£ + rZ2E + ^ I ... (42 ) 



V2 and V3 are distinguished in that they are linear combinations of 611 (KdV) and new scale-invariant 

u 2 

advecting velocities ^ or -f . V n >n differ from KdV by non-scale-invariant advecting velocities. 
They also involve higher order derivatives than in the dispersive term u xxx . 

Thus, KdV is at the center of an infinite dimensional space of advection-dispersion equations with 
sech 2 solitary wave solutions. In the space of equations, we can go out a 'distance' 1 + 8 in any of 
the directions defined by V n while retaining this solution. In the simplest (least non-linear) case V2 , 
we get the KdV-SIdV family interpolating between KdV (6 = -1) and ([3]) {8 = 1 ) 

u t + (3(1 - S)u + (S+l)^-)u x = 8 u X xx- (43) 

It may be noted that 8 and e are distinctly different parameters and we have therefore kept a separate 
notation, e is a dispersion coefficient that is exclusive to SIdV advection and 8 is a different parameter 
that allows the KdV and the SIdV(when 6=1) equations to be regarded as part of the SIdV-KdV one 



parameter family. Like KdV and SIdV, ( 43 1 preserves J u 2 dx since it can be written in conservation 
law form 

( 1 - 8)u 3 + [\+8\u 2 - 8uu x \ = 0. (44) 



As we omitted a constant of integration in (40 1, the above argument does not guarantee that every 
other KdV travelling wave is a solution of (43 1. For example, there are non-decaying KdV travelling 
waves (e.g. 2sech 2 (;t + 2t) - 1) which are not solutions of (43 1 when 5=1. Nevertheless, ^ has 



cnoidal and non-decaying sech 2 solitary waves of its own. Had we kept the constant of integration 
-cf + 3f 2 + f" = A, we would have found V = 3(1 - 8)f + (8 + 1) (— ~ j)- varying A we get 
equations that share particular asymptotic classes of travelling wave solutions with KdV. 



7 Discussion and conclusion 

In this paper we studied a scale-invariant analogue of the KdV equation, u t + 2u xx u x /u = eu xxx , 
which we named the SIdV equation. SIdV is one of the two simplest translation, scale and space- 
time parity-invariant non-linear advection-dispersion equation^] The dimensionless parameter e 
measures the strength of dispersion relative to advection. For e = 1 , SIdV shares the sech 2 solitary 
wave with KdV, as originally discovered by genetic programming. When e = ±2/3,0, the equation 
reduces to the integrable mKdV equation, to a linear dispersive wave equation and to a non-linear 
diffusion equation. Like KdV, SIdV admits similarity solutions, solitary and cnoidal travelling waves, 
and remarkably, even plane waves. SIdV may be written in conservation law form in two ways (|5][6]), 
leading to two integrals of motions. In general, we have not been able to find a Lagrangian for SIdV, 
indeed we could show that there is no Lagrangian polynomial in <p , where <p x = log u . However, we 
found a Lagrangian in a different variable at e = -2/3 and a Hamiltonian when e - 2/3 . 

For e - 1 , numerical evolution of a pair of solitary waves did not show KdV-like soliton scatter- 
ing, nor did the wave intensity tend towards equipartition among all Fourier modes, even after three 
dozen collisions. The wave intensity appeared to circulate among the modes present in the initial 
state. Effectively, a finite number of degrees of freedom appeared to take part in the dynamics. This 

12 The other one involves the reciprocal of SIdV's advecting velocity: u, + —u x = ±u xxx . 
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was also manifested in the boundedness of the gradient energy Q in our simulations. So despite 
the apparent absence of soliton scattering and presence of only two known conserved quantities, Q 
seems to display 'Birkhoff recurrence' like the Fermi-Pasta-Ulam or KdV systems. Thus, it may 
provide a counter example to the idea that integrability is necessary for recurrence. This was sug- 
gested in [SI, where the number of effective degrees of freedom of a conservative nonlinear wave 
equation was identified as the possible origin of a Birkhoff recurrence, a concept generalizing the 
well-known Poincare recurrence of finite-dimensional Hamiltonian systems. The above-mentioned 
SIdV-like equation with reciprocal advecting velocity V = u/u xx is also interesting from this view- 
point, as it admits the conserved quantity J u 2 x dx. It follows by Poincare's inequality that the mean 
square mode number Q is bounded, and we expect it too to display recurrence ! 



There are of course several other interesting directions such as perturbation theory for ( 43 1 around 
the KdV limit, and for SIdV around the solvable e = ±2/3 limits. Additional conserved quantities 
and solutions on bounded domains with homogeneous boundary conditions or on a half-line are also 
of interest. More generally, we would like to study structural issues like regularity and stability. We 
are also examining multi-dimensional analogues of SIdV as well as a complex version (i\j/ t + ^\f/ x = 
iftxxx) which has a locally conserved 'probability' density \i//\ 2 . 

Motivated by the computer-aided discovery of Q we searched for other advection-dispersion 
equations u t + Vu x = Su xxx sharing the KdV solitary wave. We found that there is an infinite sequence 



of 1 -parameter families of advective velocities V n (u;6) (41 42 1 that generalize the KdV equation 
while retaining its sech 2 solitary wave solution. They involve a linear combination of KdV-like 
V oc u and new rational advective velocities. The first two of these have the feature of being scale- 
invariant: V oc u xx /u (SIdV) and u x /u 2 . Among all these equations sharing KdV solitary waves, the 
KdV-SIdV family (43) appears to be unique in not involving higher order derivatives and admitting a 
positive definite conserved density (w 2 ). It would be interesting to study this generalization in greater 
detail. 
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A Computer-aided discovery of SIdV equation 

The basic idea behind Genetic Programming (GP) is to simulate a stochastic process by which genetic 
traits evolve in offspring, through a random combination of the genes of the parents. Following the 
seminal work by Koza [15], the GP framework provides a very useful stochastic engine to discover 
various solution regimes in a complex search terrain of a given problem. It is known that an evolu- 
tionary method is especially useful when direct methods are not available. In order to set up a GP 
engine, a non-linear chromosome structure representing a candidate solution is set up that can poten- 
tially grow to a true solution by successive applications of GP operators of selection, copy, crossover 
and mutation. The quality of a given chromosome is defined and scaled down typically to a fitness 
range [0, 1] with fitness 1 signifying a true solution. Stochastically generated chromosomes fill an 
initial pool that is evolved through successive generations in which potentially strong candidate chro- 
mosomes are selected based on their fitness values. They undergo possible refinements through GP 
operators and hopefully march towards a solution with fitness=l. In using GP to deduce a PDE like 
the KdV (in symbolic form) from an analytic travelling wave solution, one begins by considering a 
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general expression for a third order ODE, 

f"(g) = C(g, f(g), f(g), /"(£)) (45) 

where f(g) is a function of the travelling wave phase £ = x — ct . For example, a chromosome during 
a GP iteration could be C = 1.1/" + 2f'(f - 3/) + The fitness parameter is then estimated at 
each stage by examining the mean-squared difference between the 'chromosomal value' of /"' and 
its value at the given analytic solution. GP follows a 'fitness driven evolution path' by minimizing the 
error in admitting the given function as a solution. We carried out a number of GP experiments for 
the sech 2 KdV solitary wave §2^. Due to a stochastic search procedure adopted by GP, it was found to 
be too slow. We improved and accelerated it by introducing a sniffer technique Ifl2ll that carried out a 
local search at regular intervals to enhance the minimization procedure. Our improved GP approach 
was quite successful in inferring PDEs. Starting with the KdV solitary wave and its derivatives, the 
method not only reproduces the KdV equation, but also gives the SIdV equation ([3]). 
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